home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / language / embedded / mcu11 / ffthc11.arc / FFTHC11.ASM
Assembly Source File  |  1987-01-29  |  13KB  |  377 lines
  1. * fft.c11 - fast fourier transform for the MC68HC11
  2. *        written by:
  3. *
  4. *        Ron Williams
  5. *        Department of Chemistry
  6. *        Ohio University
  7. *        Athens, OH 45701
  8. *
  9. * This is a modification of the 6800 FFT presented by:
  10. *      Richard Lord
  11. *      Byte Magazine, pp. 108-119
  12. *      February 1979
  13. *
  14. * My version is written in ROMable code for the HC11. 
  15. * It uses a sine look-up table for speed and can only 
  16. * transform 256 8-bit data points.  The program 
  17. * assumes that the address of the real data is pushed 
  18. * on the stack prior to the call and that a 256 byte 
  19. * imaginary buffer is at data+256 therefore you must 
  20. * declare a 512 byte data array in the calling routine 
  21. * and load the lower 256 bytes with data.  The FFT 
  22. * will zero out the imaginary portion.  Also note that 
  23. * the FFT uses memory in the stack RAM for its dynamic 
  24. * variables and the FFT returns a value on the stack 
  25. * which contains the number of times the data was 
  26. * divided by 2 during transform. 
  27. * 
  28. * As mentioned in Lord's article, "power" spectra can 
  29. * be computed using the sum of absolute values routine 
  30. * included at the end of the FFT. Simply change the 
  31. * beq done at the end of the FFT to beq smsq. 
  32. * 
  33. * note - this copy has been modified to use $DD00 for 
  34. * data because this is easiest with BUFFALO
  35. * I have timed this transform on some test data.  The 
  36. * results are an impressive 350 milliseconds per 
  37. * transform including the "power" spectra computation. 
  38. * 
  39. * Please note that the origin of this code may require 
  40. * adjustment for your specific memory map. 
  41. * 
  42. * Please let me know of any bugs you find.
  43.  
  44. return  equ     0
  45. real    equ     2
  46. celnm   equ     4
  47. celct   equ     5
  48. pairnm  equ     6
  49. celdis  equ     7
  50. delta   equ     8
  51. sclfct  equ     9
  52. cosa    equ     $0A
  53. sina    equ     $0B
  54. sinpt   equ     $0C
  55. real1   equ     $0E
  56. real2   equ     $10
  57. treal   equ     $12
  58. timag   equ     $13
  59. tmp     equ     $14
  60. tmp2    equ     $15
  61. data    equ     $DD00
  62.  
  63.         org     $C000
  64.  
  65.         tsx             top of stack for frame pointer
  66.         xgdx                to be placed in X
  67.         subd    #$18    subtract offset to make room
  68.         xgdx            X now has frame pointer
  69.         puly            get return address
  70.         sty     return,X save it
  71. *        puly            get data address
  72. *        sty     real,X  save it
  73.         ldy    #data
  74.         sty    real,X
  75.         clr     sclfct,X zero scale factor
  76.         iny             inc y for imag data
  77.         clrb
  78. zero    clr     $FF,Y   note special place of imag
  79.         iny              256 above data
  80.         decb
  81.         bne     zero
  82. *
  83. * must do bit sorting before transforming
  84. *
  85.         ldab    #$FE    setup start for bit reversal
  86. revbit  ldaa    #08     get # of bits to reverse
  87.         pshb            save address offset
  88. rev1    rorb            rotate b right - bit to carry
  89.         rol     tmp,X   rotate left - carry bit in
  90.         deca            decrement counter
  91.         bne     rev1    go back if not done
  92.         pulb            get unshifted address
  93.         pshb            save copy
  94.         cmpb    tmp,X   check to see if already done
  95.         bcs     noswap  if so don't swap bytes
  96. swap    ldy     real,X  get data address
  97.         aby             add to base address
  98.         ldaa    0,Y     get value
  99.         pshy            store away
  100.         ldy     real,X  get base again
  101.         ldab    tmp,X   get shifted address
  102.         aby             add to base
  103.         ldab    0,Y     get second member
  104.         staa    0,Y     put away first member
  105.         puly            get first address
  106.         stab    0,Y     put second member in first slot
  107. noswap  pulb            get current address back
  108.         decb            decrement it
  109.         bne     revbit  do next if not done
  110. *
  111. * special case of first pass of FFT
  112. *
  113.         jsr     scale
  114.         ldy     real,X  set up data pointer
  115.         ldaa    #128    get number of cells
  116.         staa    tmp,X   store in temp
  117. fpss    ldaa    0,Y     get RM
  118.         ldab    1,Y     get RN
  119.         psha            make copy
  120.         aba             RM'=RM+RN
  121.         staa    0,Y     save back in data array
  122.         pula            get RM again
  123.         sba             RN'=RM-RN
  124.         staa    1,Y     put away
  125.         iny             point to next pair
  126.         iny
  127.         dec     tmp,X   decrement # cells
  128.         bne     fpss    go back if not done
  129. *
  130. * now the FFT proper for passes 2 thru N
  131. *
  132. four    ldaa    #64     # of cells is now 64
  133.         staa    celnm,X store
  134.         staa    delta,X so is delta
  135.         ldaa    #02     number of pairs is 2
  136.         staa    pairnm,X
  137.         staa    celdis,X so is distance between
  138. npass   jsr     scale   check for over-range
  139.         ldaa    celnm,X get current cell #
  140.         staa    celct,X store at cell counter
  141.         ldy     real,X
  142.         sty     real1,X get copy of data
  143. ncell   ldy     #sintab get address of sines
  144.         sty     sinpt,X save copy
  145.         ldaa    pairnm,X get current pairnm
  146. np1     psha            save pair counter
  147.         ldaa    0,Y     get cosine
  148.         ldab    64,Y    get sine
  149.         staa    cosa,X  save copy
  150.         stab    sina,X  ditto
  151.         ldy     real1,X point to top of data
  152.         ldab    celdis,X get current offset
  153.         aby             add to Y for current 
  154.         sty     real2,X copy it
  155.         ldaa    0,Y     get data point rn
  156.         psha            copy it
  157.         ldab    cosa,X  get cosine
  158.         jsr     smul    rn*cos(a)
  159.         staa    treal,X
  160.         pula            get copy of rn
  161.         ldab    sina,X  get sin(a)
  162.         jsr     smul    rn*sin(a)
  163.         staa    timag,X store imaginary tmp
  164.         iny
  165.         ldaa    $FF,Y   get imaginary data
  166.         psha            save it
  167.         ldab    sina,X  get sin(a)
  168.         jsr     smul    in*sin(a)
  169.         adda    treal,X
  170.         staa    treal,X  tr=rn*cos + in*sin
  171.         pula            get data back
  172.         ldab    cosa,X  get cosine
  173.         jsr     smul    in*cos(a)
  174.         suba    timag,X  ti=in*cos-rn*sin
  175.         staa    timag,X
  176.         ldy     real1,X
  177.         ldaa    00,Y    get rm 
  178.         tab             save a copy
  179.         adda    treal,X rm'=rm+tr
  180.         staa    00,Y    store new rm
  181.         subb    treal,X rn'=rm-tr
  182.         ldy     real2,X
  183.         stab    00,Y    store new rn
  184.         ldy     real1,X
  185.         iny
  186.         sty     real1,X save real1 for nxt
  187.         ldaa    $FF,Y   get im
  188.         tab             save copy
  189.         adda    timag,X im'=im+ti
  190.         staa    $FF,Y   put back in array
  191.         ldy     real2,X
  192.         iny
  193.         subb    timag,X in'=im-ti
  194.         stab    $FF,Y   put back in array
  195.         ldy     sinpt,X
  196.         ldab    delta,X increment sine pntr
  197.         aby
  198.         sty     sinpt,X save away
  199.         pula
  200.         deca            dec pair counter
  201.         bne     np1
  202. ar1     ldy     real1,X
  203.         ldab    celdis,X
  204.         aby
  205.         sty     real1,X
  206.         dec     celct,X
  207.         beq     ar3
  208.         jmp     ncell
  209. ar3     lsr     celnm,X  half cells
  210.         beq     smsq    done when all cells
  211.         asl     pairnm,X double pairs
  212.         asl     celdis,X twice as far apart
  213.         lsr     delta,X  delta is half
  214.         jmp     npass    one more time!
  215. done    ldaa    sclfct,X get scale factor
  216.         psha             save on stack
  217.         ldy     return,X
  218.         pshy
  219.         rts
  220. *
  221. * sum of absolute values instead of sum of squares
  222. *
  223. smsq    ldy     real,X   compute sum of "sqrs"
  224.         clra             clear byte counter
  225. sum     psha             save on stack
  226.         ldaa    0,Y      get real data point
  227.         bpl     sm1      force positive
  228.         nega
  229.         bvc     sm1      watch for $80
  230.         clra               which is really 0
  231. sm1     iny              get imaginary data
  232.         ldab    $FF,Y
  233.         bpl     sm2      force positive
  234.         negb
  235.         bvc     sm2      watch for $80 again
  236.         clrb
  237. sm2     dey              correct data pointer
  238.         aba              compute sum
  239.         staa    0,Y      save back in real
  240.         iny              inc Y for next round
  241.         pula             get byte counter
  242.         deca             done when zero
  243.         bne     sum
  244.         bra     done     lets get out of here
  245. *
  246. * subroutine for catching overscaled data
  247. *
  248. scale   ldy     real,X   start at top of data
  249.         ldab    #$FF
  250.         aby              top of data
  251.         aby              top or imag
  252.         iny              need two more
  253.         iny
  254.         ldaa    #$C0     -64
  255.         ldab    #$40     +64
  256. top     cmpa    0,Y      check for minimum
  257.         blo     nxt      if more negative fix
  258.         cmpb    0,Y      check for too big
  259.         bcs     scl      go fix it
  260. nxt     dey              bump pointer
  261.         cpy  real,X      done when both
  262.         bne     top      imag and data done
  263.         rts
  264. scl     inc     sclfct,X keep track of scale
  265.         ldy     real,X   set up pointer
  266.         ldab    #$FF
  267.         aby
  268.         aby
  269.         iny
  270.         iny
  271. scl1    ldaa    0,Y      get data
  272.         adda    #$80     make positive
  273.         lsra             divide by two
  274.         suba    #$40     put back
  275.         staa    0,Y      store away
  276.         dey              bump pointer
  277.         cpy     real,X   done when both
  278.         bne     scl1     imag and data done
  279.         rts
  280. *
  281. * the HC11 multiply must be modified to handle
  282. * negative data
  283. *
  284. smul    staa    tmp,X   copy multiplier
  285.         stab    tmp2,X  ditto multiplicand
  286.         tsta            check sign of multiplier
  287.         bpl     sk1     skip negation
  288.         nega
  289.         bvs     sko     check for $80
  290.         beq     sko     check for zero
  291. sk1     tstb            check multiplier sign
  292.         bpl     sk2
  293.         negb
  294.         bvs     sko     check for $80
  295.         beq     sko
  296. sk2     mul             do multiplication
  297.         adca    #0      8 bit conversion
  298.         asla            and correct for sine
  299.         ldab    tmp2,X  get original multiplicand
  300.         eorb    tmp,X   check for result
  301.         bpl     out
  302.         nega            result is negative
  303. out     rts
  304. sko     clra            return zero to main
  305.         rts
  306. *
  307. * now for the sine look up table
  308. *
  309. sintab            
  310.  fcb  127, 127, 127, 127, 126, 126, 126, 125, 125, 124
  311.  fcb  123, 122, 122, 121, 120, 118, 117, 116, 115, 113
  312.  fcb  112, 111, 109, 107, 106, 104, 102, 100,  98,  96
  313.  fcb   94,  92,  90,  88,  85,  83,  81,  78,  76,  73
  314.  fcb   71,  68,  65,  63,  60,  57,  54,  51,  49,  46
  315.  fcb   43,  40,  37,  34,  31,  28,  25,  22,  19,  16
  316.  fcb   12,   9,   6,   3,   0,  -3,  -6,  -9, -12, -16
  317.  fcb  -19, -22, -25, -28, -31, -34, -37, -40, -43, -46
  318.  fcb  -49, -51, -54, -57, -60, -63, -65, -68, -71, -73
  319.  fcb  -76, -78, -81, -83, -85, -88, -90, -92, -94, -96
  320.  fcb  -98,-100,-102,-104,-106,-107,-109,-111,-112,-113
  321.  fcb -115,-116,-117,-118,-120,-121,-122,-122,-123,-124
  322.  fcb -125,-125,-126,-126,-126,-127,-127,-127,-127,-127
  323.  fcb -127,-127,-126,-126,-126,-125,-125,-124,-123,-122
  324.  fcb -122,-121,-120,-118,-117,-116,-115,-113,-112,-111
  325.  fcb -109,-107,-106,-104,-102,-100, -98, -96, -94, -92
  326.  fcb  -90, -88, -85, -83, -81, -78, -76, -73, -71, -68
  327.  fcb  -65, -63, -60, -57, -54, -51, -49, -46, -43, -40
  328.  fcb  -37, -34, -31, -28, -25, -22, -19, -16, -12,  -9
  329.  fcb   -6,  -3,   0,   3,   6,   9,  12,  16,  19,  22
  330.  fcb   25,  28,  31,  34,  37,  40,  43,  46,  49,  51
  331.  fcb   54,  57,  60,  63,  65,  68,  71,  73,  76,  78
  332.  fcb   81,  83,  85,  88,  90,  92,  94,  96,  98, 100
  333.  fcb  102, 104, 106, 107, 109, 111, 112, 113, 115, 116
  334.  fcb  117, 118, 120, 121, 122, 122, 123, 124, 125, 125
  335.  fcb  126, 126, 126, 127, 127, 127
  336. *
  337. * Some routines for fast A2D and for sending data to 
  338. * host
  339. *   these have not been debugged
  340. *
  341. chan4   equ     $04
  342. atdctr  equ     $1030
  343. portc   equ     $1003
  344. ccontr  equ     $1002
  345. lportc  equ     $1005
  346. strt    equ     $02
  347.         ldy     #data
  348.         ldaa    #chan4
  349.         ldab    #strt
  350. wait:   bitb    portc
  351.         beq     wait
  352. cklp:   ldab    ccontr
  353.         bpl     cklp
  354.         staa    atdctr
  355.         ldab    lportc
  356.         iny
  357.         cpy     #data+256
  358.         beq     adone
  359.         nop
  360.         nop
  361.         nop
  362.         nop
  363.         ldab    atdctr+1
  364.         stab    $00,Y
  365.         bra     cklp
  366. adone:  rts
  367.  
  368. output  equ    #$E3B3
  369. sendat  ldx    #data
  370.         clrb
  371. s1      jsr    output
  372.         decb
  373.         bne    s1
  374.         rts
  375.  
  376.  
  377. ut  equ    #$E3B